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More  efficient  water  management  techniques  are  required  to  decrease  the  cost  of  polymer  electrolyte 
fuel  cell  (PEFC)  systems  while  maintaining  robust  performance.  In  this  study,  we  use  nonlinear  statistical 
analysis  of  experimental  data  to  characterize  PEFC  dynamics  under  conditions  where  water  accumulation 
in  the  cathode  air-delivery  microchannels  causes  decreases  in  performance  accompanied  by  chaotic 
fluctuations.  Using  experimental  PEFC  voltage  signals,  we  estimate  chaotic  invariants  indicative  of  the 
degrees  of  freedom  of  the  dynamics  (the  correlation  dimension)  and  the  instability  of  the  dynamics  (the 
Kolmogorov  entropy).  We  find  that  these  invariants  decrease  with  increasing  gas  flow  commensurate 
with  greater  fuel  cell  current  and  air  stoichiometric  ratio,  and  that  they  are  indicative  of  the  channel  two- 
phase  flow  regime.  We  correlate  the  Lyapunov  exponents  of  the  one-dimensional  voltage  return  map  and 
the  Hurst  exponents  of  the  voltage  time  series  with  the  chaotic  invariants  for  use  in  future  PEFC  water 
management  and  control  strategies.  In  addition,  we  examine  the  relationships  between  the  invariants 
estimated  from  the  voltage  signal  and  the  two-phase  friction  multiplier  calculated  from  measured 
cathode  pressure  drop  in  order  to  distinguish  the  distinct  dynamics  of  two-phase  channel  flow. 

©  2014  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Polymer  electrolyte  fuel  cells  (PEFCs)  are  electrochemical  en¬ 
gines  that  convert  hydrogen  and  oxygen  to  electricity,  with 
byproducts  of  water  and  heat.  The  ability  of  sustainably-  and 
renewably-produced  hydrogen  to  provide  low  emission  electricity 
makes  it  an  attractive  alternative  to  the  combustion  of  limited- 
supply,  carbon-based  fuels  like  petroleum  and  coal.  Because 
PEFCs  offer  high  power  density  and  quick  startup  owing  to  their  sub 
100  °C  operating  temperature,  they  remain  promising  for  use  in 
transportation  applications. 


*  Corresponding  author. 

E-mail  addresses:  mburkhol@andrew.cmu.edu  (M.B.  Burkholder),  nicholas. 
siefert@netl.doe.gov  (N.S.  Siefert),  litster@andrew.cmu.edu  (S.  Litster). 

http://dx.doi.org/10.1016/jjpowsour.2014.04.156 

0378-7753 /©  2014  Elsevier  B.V.  All  rights  reserved. 


While  the  product  water  exhaust  makes  PEFCs  attractive  as 
clean  engines,  it  also  poses  water  management  problems  that 
reduce  system  efficiency  and  hinder  robustness  of  operation  [1  .  In 
the  PEFC  cathode,  water  is  produced  according  to  the  oxygen 
reduction  reaction  in  which  oxygen,  usually  provided  by  air  flow, 
combines  with  hydrogen  ions  and  electrons.  A  certain  amount  of 
water  is  desirable  to  hydrate  the  polymer  electrolyte  membrane, 
thereby  decreasing  its  ionic  resistance  [2  .  However,  too  much 
water  can  increase  the  resistance  to  oxygen  diffusion,  leading  to 
oxygen  starvation  and  large  mass  transport  overpotentials. 

The  ability  to  effectively  remove  surplus  water  from  the  cathode 
with  a  low  parasitic  cost  is  confounded  by  the  air  delivery  mech¬ 
anism  and  two-phase  flow  effects  under  certain  operating  condi¬ 
tions.  State  of  the  art  air  delivery  flow  fields  are  composed  of  arrays 
of  parallel  microchannels  which  deliver  oxygen  to  the  active  fuel 
cell  area  and  remove  excess  water.  To  achieve  low  parasitic  costs,  it 
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is  desirable  to  use  low  air  flow  rates  in  channels  with  low  hydraulic 
resistance.  However,  the  air  momentum  that  makes  positive  con¬ 
tributions  to  the  hydraulic  resistance  also  removes  liquid  water 
through  nonlinear  momentum  transfer.  Hence,  low  air  flow  rates 
allow  water  to  accumulate  in  the  microchannels,  thereby  reducing 
the  effective  cathode  active  area  and  the  fuel  cell  power  output. 

For  high  current  densities,  neutron  radiography  visualization 
has  shown  that  the  water  accumulation  in  the  cathode  micro¬ 
channels  is  relatively  insensitive  to  the  air  stoichiometric  ratio  £  [3], 
which  is  the  ratio  of  the  air  flow  rate  delivered  to  that  required  for 
the  reaction.  At  low  current  densities,  however,  the  water  accu¬ 
mulation  and  PEFC  performance  become  more  sensitive  to  f 
because  the  two-phase  flow  regime  in  the  microchannels  tends 
towards  slug/plug  flow  when  the  channel  superficial  air  and  water 
velocities  (ug  and  u/,  respectively)  are  low  [4,5  .  In  arrays  of  parallel 
channels,  the  complete  blockage  of  individual  channels  by  water 
plugs  allows  air  flow  to  divert  from  plugged  channels  to  open 
channels,  following  the  path  of  least  resistance  [6  .  The  diverted 
flow  leaves  flooded  channels  plugged  until  enough  pressure  builds 
up  to  clear  them,  resulting  in  voltage  instability  and  fluctuations, 
and  power  loss  [7,8]  (see  Fig.  1 ).  While  the  channel  plurality  [8]  in  a 
single  lab-scale  fuel  cell  may  prove  manageable,  the  large  numbers 
of  parallel  channels  present  in  fuel  cell  stacks  used  to  power  an 
automobile  (>10,000)  introduce  many  additional  degrees  of 
freedom  and  instability  into  the  system's  dynamics. 

To  avoid  detrimental  instability  from  water  accumulation, 
operation  in  the  slug  flow  regime  is  minimized  in  favor  of  annular 
film  flow  [5  ,  which  can  be  aided  through  geometric  and  material 
flow  field  design  optimization.  Different  geometric  flow  field 
designs  have  been  investigated,  and  while  interdigitated  chan¬ 
nels  have  been  shown  to  effectively  remove  liquid  water  9],  they 
produce  much  larger  pressure  drops  than  the  conventional  par¬ 
allel/serpentine  channel  flow  fields.  To  promote  stable  liquid 
water  removal  via  annular  film  or  corner  flow,  hydrophilic  flow 
field  materials  may  be  used  [10,11  as  well  as  smaller  channel 
dimensions  that  increase  the  gas  velocity.  Even  with  these  flow 
field  design  optimizations,  effective  water  removal  is  driven  by 
supplying  excessive  air  flow  rates  at  multiple  times  the 


Fig.  1.  Voltage  signal  of  an  operating  PEFC  beginning  from  a  purged  condition.  As 
water  accumulates,  voltage  output  decreases  and  chaotic  voltage  fluctuations  are 
introduced  (around  400  s)  before  the  cathode  completely  floods  and  voltage  is  lost 
(around  700  s). 


stoichiometric  flow  rate,  resulting  in  high  parasitic  load  on  the 
fuel  cell  output,  especially  during  the  low-power  idling  and 
startup  conditions. 

One  class  of  PEFC  water  management  techniques  uses  fault 
detection.  Output  indicative  of  fuel  cell  health  is  monitored,  and 
corrective  water  purging  action  is  taken  when  the  output  indicates 
that  the  cathode  is  flooded.  Many  different  diagnostics  have  been 
proposed  for  fault  detection.  One  such  linear  fault  detection 
scheme  monitors  the  levels  and  oscillation  frequencies  of  fuel  cell 
voltage  and  resistance,  stack  voltage,  pressure  drop,  and  load  12]. 
When  any  of  these  signals  breaches  a  predefined  threshold,  the  air 
flow  rate  is  pulsed  to  drive  out  liquid  water.  The  amplitude 
response  of  voltage  to  small  output  pressure  oscillations,  indicated 
by  the  voltage  variance,  has  also  been  proposed  as  a  diagnostic  for 
cathode  flooding  [13  .  Additionally,  electrochemical  impedance 
spectroscopy  has  been  employed  to  indicate  cathode  health  14]. 
Fuel  cell  dehydration  shows  a  strong  impedance  response  across  a 
broad  frequency  range,  while  cathode  flooding  shows  a  much 
weaker  impedance  response  across  a  limited  frequency  range.  Both 
of  the  last  two  techniques  require  additional  hardware  to  periodi¬ 
cally  excite  the  PEFC  system,  with  the  latter  also  requiring  fre¬ 
quency  response  analysis  electronics. 

To  improve  upon  such  fault  detection  techniques,  the  signal 
being  monitored  and  the  threshold  values  can  be  optimized  in  a 
way  that  is  consistent  with  the  nonlinearities  of  PEFC  water  accu¬ 
mulation  and  removal  dynamics.  Nonlinear  chaotic  behavior  has 
long  been  observed  in  chemical  reaction  kinetics,  like  the  classic 
Belousov-Zhabotinsky  reaction  15].  The  complex  PEFC  water 
balance  involves  two-phase  flow  in  parallel  microchannels  with 
electrochemical  feedback,  resulting  in  highly  nonlinear  and  chaotic 
dynamics.  While  linear  indicators  like  voltage  oscillation  frequency 
and  magnitude  may  indicate  changes  in  the  PEFC  dynamics,  they  do 
not  account  for  potential  nonlinear  differences  (e.g.  in  stability) 
between  dynamics  with  similar  linear  oscillation  characteristics.  To 
address  the  shortcomings  of  linear  analysis,  nonlinear  techniques 
have  been  used  to  analyze  similar  chaotic  systems.  For  instance, 
calculation  of  chaotic  invariants  like  Kolmogorov  entropy  has  been 
shown  to  indicate  two-phase  flow  regime  changes  16],  and 
nonlinear  statistics  have  been  shown  to  provide  early  defluidiza¬ 
tion  warning  in  a  fluidized  bed  [17  .  Thus,  instability  induced  by 
changes  in  system  parameters  with  long  dynamical  timescales, 
such  as  fuel  cell  temperature  and  pressure,  is  potentially  detectable 
using  nonlinear  statistics  before  it  is  apparent  from  linear  statistics. 
It  is  also  possible  to  stabilize  unstable  dynamics  with  minimal  effort 
using  chaos  control  techniques  without  changing  the  system  pa¬ 
rameters  [18].  Bubbling  has  been  controlled  with  chaos  control  in  a 
fluidized  bed  [17  ,  and  chaos  has  been  observed  and  controlled  in 
electrochemical  reactions  [19  . 

The  purpose  of  this  paper  is  to  provide  a  framework  for  un¬ 
derstanding  the  nonlinear,  chaotic  dynamics  of  PEFC  water  man¬ 
agement.  We  apply  the  first  nonlinear  time  series  analysis  to  PEFC 
instabilities  to  investigate  how  the  presence  of  two-phase  flow 
affects  invariant  quantities  indicative  of  the  degrees  of  freedom 
and  the  stability  of  the  dynamics.  We  also  propose  new  nonlinear 
voltage  statistics  that  indicate  the  dynamical  stability  and  are 
much  less  expensive  to  compute  than  the  traditional  chaotic 
invariants. 

2.  Method 

This  section  presents  the  nonlinear  analysis  methodology  that 
we  used  in  this  work  and  insight  into  its  physical  interpretation. 
Our  implementation  of  the  equations  and  algorithms  has  been 
validated  on  well-characterized  nonlinear  and  stochastic  systems 
like  the  Lorenz  system  [20]  and  Brownian  motion,  respectively. 
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2.2.  State  space 

Due  to  a  lack  of  analytical  tools,  common  analysis  techniques  for 
nonlinear  systems  focus  on  the  geometrical  properties  of  the  sys¬ 
tem  state-space  attractor.  When  nonlinear  differential  state  equa¬ 
tions  are  integrated,  they  yield  a  trajectory  through  the  state  space 
over  time,  where  each  orthogonal  coordinate  on  the  trajectory  is  a 
state  as  a  function  of  time.  The  region  of  state  space  to  which  the 
trajectory  is  confined  is  the  system  attractor. 

Experimentally,  the  state  space  can  be  reconstructed  from  the 
observed  time  series  of  a  single  state,  x(t)  [21  ].  The  technique,  called 
time-delay  embedding,  creates  a  time  series  of  m-dimensional 
embedded  vectors,  x(t),  where  each  coordinate  is  the  observed  time 
series  lagged  a  time  delay  r  behind  the  previous  coordinate,  Eq.  (1 ). 
Such  a  reconstructed  vector  space  preserves  the  topology  of  the 
system  state  space,  and  hence  the  geometrical  properties  of  the 
attractor,  provided  that  enough  dimensions,  in  this  case  time- 
delays,  are  used. 

x(t)  =  [x(t),x(t  -  T),x(t  -  2t),  -  (m  -  1  )t)]  (1) 

Using  the  embedded  space,  invariant  properties  of  the  state 
space,  i.e.  quantities  that  are  the  same  for  all  valid  state  space  re¬ 
constructions,  may  be  estimated  with  neighborhood  statistics. 
Neighborhood  statistics  are  used  to  determine  not  only  how  a 
single  point  on  the  trajectory  evolves  over  time,  but  also  how 
similar  states  within  a  neighborhood  of  distance  r  evolve  over  time. 
For  non-chaotic  systems,  the  evolution  of  similar  states  is  also 
similar  and  neighborhood  statistics  provide  little  additional  infor¬ 
mation.  In  chaotic  systems,  due  to  the  sensitive  dependence  on 
initial  conditions,  similar  states  exhibit  exponentially  differing  time 
evolution  and  thus  neighborhood  statistics  are  necessary  to  quan¬ 
tify  attractor  properties.  This  is  because  the  evolution  of  a  single 
initial  condition  is  not  indicative  of  the  evolution  of  other  initial 
conditions,  even  similar  ones.  Fig.  2  shows  a  three-dimensional 
section  of  the  >3D  chaotic  time-delay  embedded  attractor  for  the 
voltage  signal  in  Fig.  1. 

Confidently  estimating  state  space  invariants  requires  accurate 
neighborhood  statistics  over  a  large  range  of  neighborhood  sizes. 
For  an  experimentally  measured,  time-delay  embedded  state  space, 
finite  amounts  of  data  can  limit  statistical  accuracy  at  large  r.  In 


Fig.  2.  Time-delay  embedded  attractor  of  the  voltage  signal  in  Fig.  1  with  m  =  3.  The 
strange  attractor  in  the  chaotic  region  has  a  fractal  dimension. 


addition,  the  presence  of  dynamical  and  measurement  noise  scat¬ 
ters  the  points  around  the  attractor  manifold,  limiting  accuracy  at 
small  r.  Due  to  the  computational  intensity  of  the  estimation  al¬ 
gorithms,  it  is  often  impractical  to  use  large  amounts  of  data, 
making  noise  reduction  desirable  to  increase  the  accessible  range  of 
r. 

Low-pass  filters  commonly  used  for  noise  reduction  are  based 
on  the  assumption  that  frequency  components  from  noise  are  lin¬ 
early  additive  to  the  signal.  In  nonlinear  systems,  power  input  at 
one  frequency  is  distributed  among  many  frequencies,  violating 
this  assumption.  Thus,  we  used  a  nonlinear  noise  reduction  (NNR) 
algorithm  [22]  to  filter  the  voltage  signals.  The  noise  reduction, 
schematically  depicted  in  two  dimensions  in  Fig.  3,  is  performed  in 
the  reconstructed  state  space.  For  each  point  in  the  m-dimensional 
space  (m  =  2  in  Fig.  3),  a  hyper  ellipsoid  is  fit  to  the  neighborhood. 
Noise  is  reduced  by  projecting  the  point  towards  the  local,  linear¬ 
ized  attractor  manifold  approximated  by  the  neighborhood.  The  q 
smallest  principal  axes  of  the  hyper  ellipsoid  are  assumed  to  be 
noise  components  (q  =  1  in  Fig.  3),  and  each  point  is  projected 
toward  the  m-q- dimensional  attractor  manifold  by  removing  the 
noise  components  (m-q  =  1  in  Fig.  3).  The  procedure  is  iterated  for 
convergence. 

2.2.  Fractal  dimensions 

One  invariant  of  the  state  space  that  we  estimated  is  the 
attractor  dimension,  which  is  the  active  degrees  of  freedom  avail¬ 
able  to  the  dynamics.  The  attractor  dimension  typically  has  a  fractal 
value  in  the  presence  of  chaotic  dynamics,  and  for  this  reason 
chaotic  attractors  are  often  called  strange  attractors. 

The  topological  dimension  is  an  integer  describing  how  many 
coordinates  are  necessary  to  specify  a  unique  point  in  any  object, 
including  dynamical  attractors.  The  topological  dimension  of  an 
attractor  is  how  many  effective  states  contribute  to  the  dynamics, 
or  the  number  of  state  space  dimensions  in  which  the  system  dy¬ 
namics  are  free  to  move.  A  stable  m-dimensional  state  space  tra¬ 
jectory  will  typically  be  confined  to  an  integer-dimensional 
attractor  (of  dimension  n)  in  a  subspace  of  the  state  space  (n  <  m), 
with  a  higher  n  being  associated  with  more  complex  dynamics. 
Should  a  manifold  in  the  state  space  become  unstable,  as  through  a 


Fig.  3.  Schematic  of  the  nonlinear  noise  reduction  procedure,  wherein  each  point  is 
projected  towards  the  attractor  manifold  approximated  by  the  primary  axes  of  the 
neighborhood  best-fitting  ellipsoid.  Here,  the  noise  component  of  the  square  point  is 
depicted,  as  well  as  the  trajectory  before  and  after  the  noise  reduction  procedure. 
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bifurcation,  the  dynamics  become  chaotic  and  the  attractor  gains 
additional  degrees  of  freedom  in  the  direction  of  instability, 
resulting  in  a  fractally-dimensioned  attractor.  From  the  point  of 
view  of  the  integer  topological  dimension,  the  chaotic  attractor 
does  not  entirely  fill  another  dimension  of  state  space  and  thus  still 
has  a  dimension  of  n,  when  in  reality  the  fractal  attractor  fills  more 
space  than  it's  topological  dimension  of  n  but  less  space  than  n  +  1 
dimensions. 

For  this  reason  it  is  necessary  to  generalize  the  topological 
dimension  to  account  for  the  additional  degrees  of  freedom  of  the 
dynamics.  One  such  generalization  is  the  correlation  dimension,  D2, 
a  straightforward  estimate  of  attractor  dimension  in  a  time-delay 
embedded  state  space  [23].  Estimation  of  D2  is  achieved  by 
computing  correlation  sums  C(r)  from  Eq.  (2),  which  represent  how 
the  average  attractor  density  (or  number  of  neighbors)  changes  as  a 
function  of  the  neighborhood  radius  r.  The  power  law  scaling  of  C(r) 
with  r  then  yields  the  correlation  dimension:  C  °c  r°2 ,  which  gives  an 
integer  value  for  attractors  with  integer  dimensions  and  a  decimal 
value  for  attractors  with  fractal  dimensions.  In  Eq.  (2),  8  is  the 
Heaviside  step  function,  N  is  the  total  number  of  points,  and  since 
we  don't  want  to  count  points  that  are  in  the  neighborhood  due  to 
temporal  correlations,  T  is  the  length  of  the  temporal  correlation.  By 
averaging  the  number  of  points  within  the  neighborhood  of  size  r, 
we  then  have  the  correlation  sums  as  a  function  of  r,  and  also  as  a 
function  of  m  through  the  embedding  procedure. 
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Due  to  the  finite  amount  of  data  and  sampling  frequency,  the 
power  law  scaling  necessary  to  make  an  estimate  of  D2  occurs  over 
a  limited  range  of  r,  called  the  scaling  range.  A  dimension  estimate 
should  only  be  made  when  there  is  clear  linearity  in  the  log-log 
plot  of  C(r),  which  is  apparent  from  a  plateau  in  the  local  slopes. 
Because  noise  can  also  produce  a  plateau  in  C(r)  at  the  embedding 
dimension  m,  a  plateau  for  a  confident  D2  estimate  should  be 
witnessed  for  all  correlation  sums  with  embedding  dimensions 
high  enough  to  accommodate  the  attractor. 

Fig.  4  shows  an  example  of  estimating  D2  from  a  PEFC  voltage 
signal.  The  power  law  scaling  of  the  correlation  sums  yields  a  D2 
estimate  of  approximately  10  for  a  current  density  of  0.1  A  cm”2  and 
?  =  1.5,  indicating  the  presence  of  high-dimensional  chaos.  In  our 
analysis,  we  normalized  the  delay  vectors  by  the  standard  deviation 
of  the  signal  and  so  r  is  the  unitless  number  of  standard  deviations. 


2.3.  Information  entropy 


r  [-] 

Fig.  4.  a)  Voltage  time  series  at  0.1  A  cm-2,  £  =  1.5;  b)  correlation  sums  computed  from 
the  normalized  voltage  time  series  for  m  =  18,  20,  22,  24,  26,  and  28;  c)  local  slopes  of 
the  correlation  sums  yielding  a  D2  estimate  of  approximately  10;  d)  Kolmogorov  en¬ 
tropy  estimates  from  the  correlation  sums  yielding  a  I<2  estimate  of  about  0.6  nats  s-1. 
Such  estimates  may  only  be  made  from  convergence  in  the  power  law  scaling  range 
indicated  by  a  plateau  in  the  correlation  sum  slopes. 


Another  invariant  that  we  estimate  is  the  Kolmogorov  entropy. 
It  is  an  entropy  in  the  information  theory  sense  in  the  fact  that  it  is  a 
measure  of  the  amount  of  information  about  prior  states  needed  to 
specify  the  current  state.  In  a  stable,  non-chaotic  system  where  the 
states  evolve  in  a  predictable  manner,  no  Kolmogorov  entropy  is 
present.  For  random  systems  where  no  amount  of  prior  informa¬ 
tion  can  be  used  to  predict  the  next  outcome,  there  is  infinite 
Kolmogorov  entropy.  For  chaotic  systems  with  sensitive  depen¬ 
dence  on  initial  conditions  due  to  the  exponential  divergence  of 
initially  similar  states,  there  is  a  positive,  finite  Kolmogorov 
entropy. 

The  Kolmogorov  entropy  I<  is  computed  by  dividing  the  state 
space  into  hypercubes  and  considering  the  joint  probabilities  with 
which  the  trajectory  visits  the  hypercubes.  With  an  experimental 
time  series,  it  is  much  easier  to  calculate  the  order-2  Renyi  entropy 
K 2,  which  is  positive  for  chaotic  systems  and  is  upper  bounded  by  I< 


[24].  The  computation,  described  by  Eq.  (3),  can  be  made  directly 
from  the  correlation  sums  for  varying  embedding  dimensions  m, 
where  n  is  the  integer  distance  between  consecutive  values  of  m. 


lim 

m— >00 


lim  — 

r->0  Hr 


log (  C(r’m>  ) 

6  \C(r,  m  +  n)J 


Similar  to  D2  estimation,  a  confident  I< 2  estimate  should  only  be 
made  over  the  scaling  range  when  power  law  scaling  is  apparent  in 
the  correlation  sums.  While  the  power  law  scaling  itself  is  not  a 
sufficient  condition  to  indicate  chaos,  a  positive  K 2  value  such  as 
that  shown  in  Fig.  4  is  sufficient  [24  . 

The  Kolmogorov  entropy  is  also  the  sum  of  the  positive  Lyapu¬ 
nov  exponents  (LEs,  A).  The  Lyapunov  exponent  spectrum  is  the 
exponential  rates  of  convergence  (negative)  or  divergence 
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(positive)  over  time  of  a  small  volume  of  initial  conditions  on  the 
attractor.  The  spectrum  has  as  many  exponents  as  the  attractor 
dimension,  but  usually  only  the  positive  exponents  which  charac¬ 
terize  the  level  of  instability  in  the  system  are  of  interest.  Chaotic 
systems  have  at  least  one  positive  Lyapunov  exponent,  which  gives 
rise  to  the  sensitive  dependence  on  initial  conditions  through 
exponential  divergence  of  initially  similar  states. 

The  units  of  entropy  depend  on  the  base  of  the  logarithm  used  in 
Eq.  (3).  Herein,  we  have  used  the  natural  logarithm,  making  our 
entropy  estimates  in  units  of  nats  s-1,  where  a  nat  is  the  natural  unit 
of  entropy  or  information.  Comparatively,  we  have  estimated  the 
LEs  from  a  map  with  a  base-2  logarithm  (see  next  section),  yielding 
units  of  bits  iteration1.  One  nat  is  approximately  1.44  bits. 

2.4.  Maps 

As  mentioned  earlier,  when  performing  certain  neighborhood 
statistics  to  determine  invariant  properties  of  the  attractor,  we  do 
not  consider  points  that  are  within  the  neighborhood  because  of 
temporal  correlations.  These  are  points  that  lie  on  the  same  tra¬ 
jectory  pass  through  the  neighborhood  rather  than  on  a  prior  or 
returning  pass.  In  Eq.  (2),  the  temporal  correlation  time  is  given  by 
TAt,  where  At  is  the  sampling  interval.  We  are  not  interested  in 
points  within  the  temporal  correlation  time  because  they  contain 
redundant  information,  i.e.  invariant  attractor  quantities  are 
exhibited  in  how  fast  and  how  close  trajectories  return  to  the 
original  neighborhood  after  they  have  left  it. 

A  common  technique  used  in  nonlinear  analysis  to  eliminate 
redundant  data  is  to  reduce  the  order  of  the  state  space  by  con¬ 
verting  the  time  series  ‘flow’  to  a  map.  A  map  is  a  rule  of  the  type 
xn+i  =  /(xn),  and  the  goal  is  to  choose  the  map  such  that  the  dis¬ 
carded  data  between  xn  and  xn+\  is  redundant.  A  popular  map 
choice  is  the  Poincare  map,  which  is  created  by  choosing  the  xn  as 
the  intersections  of  the  state  space  trajectory  with  a  section  in  the 
state  space.  The  choice  of  section,  called  the  Poincare  section,  thus 
defines  the  mapping,  and  the  dynamics  of  the  Poincare  map  are 
confined  to  the  Poincare  section,  a  lower-dimensional  manifold 
than  the  attractor. 

One  choice  of  Poincare  section  is  to  use  the  consecutive  maxima 
or  minima  in  the  time  series  x(t)  as  the  xn.  This  one-dimensional 
map  is  called  a  return  map.  Using  the  return  map,  parameter 
choices  like  the  time  delay  r  and  temporal  correlation  length  T  are 
eliminated.  Furthermore,  estimation  of  the  maximal  Lyapunov 
exponent  from  the  divergence  of  initially  close  peaks,  i.e.  peaks  of 
similar  magnitude,  becomes  much  more  straightforward.  For  two 
initially  similar  peaks  separated  by  a  distance  of  do,  the  local  Lya¬ 
punov  exponent  is  given  by  A  =  log2di/d0  where  d\  is  the  distance 
between  the  peaks  after  one  iteration  on  the  map  (see  Fig.  5).  The 
average  of  the  local  Lyapunov  exponents  of  many  initially  similar 
peaks  can  then  be  used  to  estimate  the  maximal  Lyapunov  exponent 
of  the  map.  Later  in  this  paper,  we  demonstrate  the  correlation 
between  the  Lyapunov  exponent  estimate  of  the  voltage  return  map 
and  the  dynamical  entropy  of  the  time-delay  embedded  state  space. 

2.5.  Diffusion  analysis 

Another  way  to  quantify  fractal  structure  in  a  signal  is  through 
computation  of  the  Hurst  exponent  [25  ,  which  was  originally 
formulated  to  study  the  long-term  storage  capacity  of  reservoirs  fed 
by  unpredictable  sources  like  rainfall  and  has  since  been  used  to 
detect  instability  in  financial  markets  prior  to  crashes  [26  ,  among 
other  applications.  For  diffusion-like  processes,  the  Hurst  exponent 
H  describes  the  power  law  scaling  between  the  distance  over  which 
the  signal  x  diffuses  and  time  t,  i.e.  (Ax2)  °c  A t2H,  where  the  angle 
brackets  denote  an  average.  For  the  special  case  where  H  =  0.5,  the 


Time  [s] 

Fig.  5.  Voltage  return  map  (dots)  for  20  s  of  the  0.1  A  cm  2,  £  =  1.5  voltage  signal.  The 
local  Lyapunov  exponent  is  estimated  from  the  divergence  of  initially  similar  peaks. 

diffusion  scales  with  the  square  root  of  time  and  is  an  uncorrelated 
random  walk  (i.e.  Brownian  motion).  A  value  of  H  different  from  0.5 
indicates  some  memory  in  the  signal,  with  H  >  0.5  indicating 
positive  autocorrelation  and  H  <  0.5  indicating  negative  autocor¬ 
relation.  A  value  of  H  =  1  is  yielded  by  a  smooth  signal  typical  of 
ballistic  motion  [27  ,  while  H  =  0  is  given  by  signals  where  adjacent 
values  lie  on  opposite  sides  of  the  mean. 

Estimation  of  the  Hurst  exponent  from  a  stationary  signal  may 
be  achieved  using  diffusion  analysis.  For  this  purpose,  a  bounded 
signal  x(t),  such  as  PEFC  voltage,  is  interpreted  as  diffusion  in¬ 
crements  to  extend  the  scaling.  Thus  the  signal  is  first  integrated  to 

yield  a  diffusion-like  quantity:  y(t)  =  [  x(t')dt'.  As  shown  in  Eq. 

Jo 

(4),  the  root  mean  squared  distance  over  which  the  integrated 
signal  diffuses  (F)  is  computed  for  varying  time  windows  of  length 
r,  where  the  angle  brackets  indicate  an  average  over  consecutive 
windows.  Power  law  scaling  between  F  and  r  gives  the  Hurst 
exponent  estimate:  F  °c  rH. 

F(r )  =  \]((y{t  +  T)  -y(t))2)  (4) 

These  analysis  techniques  for  time  series  are  nonlinear  because 
they  make  no  assumptions  about  linearity  in  the  physics  that 
produce  the  time  series.  When  we  reconstruct  an  embedded  state 
space  from  a  PEFC  voltage  signal,  we  preserve  the  invariant 
nonlinear  dynamical  properties  of  the  entire  PEFC  system, 
including  the  two-phase  flow  dynamics  of  cathode  water  accu¬ 
mulation  and  removal.  The  dimension  and  entropy  invariants,  as 
well  as  the  Lyapunov  exponents  of  the  return  map,  are  all  related  to 
the  overall  stability  of  the  dynamics,  a  property  of  concern  for 
optimal  control.  While  the  PEFC  stability  may  have  some  correla¬ 
tion  with  linear  statistics  like  voltage  mean  and  variance,  there  is  no 
guarantee  that  the  mapping  is  one-to-one.  Hence,  we  estimate 
these  invariants  for  a  PEFC  to  characterize  how  the  dynamical 
stability  changes  with  operating  conditions.  We  also  estimate  the 
Hurst  exponent  as  a  nonlinear  indication  of  the  autocorrelation,  or 
memory,  of  the  dynamics. 

3.  Experimental 

We  gathered  data  from  50  cm2  active  area  PEFC  hardware  (Fuel 
Cell  Technologies,  Inc.,  Albuquerque,  NM)  where  the  cathode  flow 


248 


M.B.  Burkholder  et  al.  /  Journal  of  Power  Sources  267  (2014)  243—254 


field  plate  was  custom  designed  and  machined  from  a  carbon- 
loaded  vinyl  ester  (BMC  940-8649,  Bulk  Molding  Compounds, 
Inc.,  Chicago,  IL)  with  a  fully  parallel  channel  flow  field  to  add  de¬ 
grees  of  freedom  and  enhance  instabilities  from  low  air  velocities 
(Fig.  6).  This  allows  better  simulation  of  full-scale  stack  dynamics 
than  a  traditional  lab-scale  serpentine  flow  field.  The  cathode  plate 
had  52  parallel  channels  covering  the  7  cm  by  7  cm  active  area,  each 
0.5  mm  wide  by  0.3  mm  deep,  corresponding  to  a  channel  plurality 
of  x  =  1  cm-2  (x*  =  1.6  x  10  3)  [8  .  Automotive  stacks  are  expected 
to  have  values  in  the  range  of  x*  ~  1.3  x  10~3  [28].  The  inlet  and 
outlet  header  dimensions  were  much  larger  than  the  channel  di¬ 
mensions  to  minimize  flow  maldistribution  from  pressure  drop 
along  the  header.  The  headers  had  pressure  taps  inserted  flush  with 
the  sidewalls  to  record  pressure  measurements  across  the  cathode 
flow  field  with  a  differential  pressure  transducer  (PX-409,  Omega 
Engineering,  Inc.,  Stamford,  CT).  The  anode  flow  field  was  a  3-pass 
serpentine  design  (see  Fig.  1  in  Ref.  [7]). 

The  membrane-electrode  assembly  was  an  18  pm  Gore  mem¬ 
brane  with  0.4/0.4  mg  Pt  cm-2  catalyst  loading  (W.  L.  Gore  and 
Associates,  Inc,  Elkton,  MD).  The  gas  diffusion  layer  was  MRC  U105 
carbon  paper  (Mitsubishi  Rayon  Co,  Ltd,  Tokyo,  Japan)  with  a 
microporous  layer. 

Fuel  cell  testing  was  conducted  using  a  commercial  fuel  cell  test 
stand  (FCT-150S,  Bio-Logic  SAS,  Claix,  France)  with  temperature, 
humidity,  load,  and  reactant  gas  flow  control.  The  fuel  cell  and 
reactant  gas  temperatures  were  held  constant  at  55  °C  for  all  ex¬ 
periments  to  achieve  100%  relative  humidity  (RH).  The  temperature 
was  lower  than  the  typical  80  °C  to  mimic  the  severe  flooding 
observed  during  startup.  The  fuel  cell  backpressure  was  atmo¬ 
spheric.  Before  each  experiment,  residual  liquid  water  was  purged 
from  the  fuel  cell  with  1000  mlpm  reactant  flow  rates  held  for 
2  min.  The  fuel  cell  was  allowed  >10  min  to  equilibrate  before  the 
data  was  gathered,  and  the  anode  flow  rate  was  held  constant  at  a 
stoichiometric  ratio  of  2. 

Galvanostatic  voltage  data  and  differential  pressure  data  were 
gathered  for  current  densities  ranging  from  0.1  to  0.3  A  cm-2  and 
air  stoichiometric  ratios,  ?,  from  1.25  to  2  at  a  frequency  of  100  Hz 
using  a  data  acquisition  board  (>10  GQ  input  impedance)  con¬ 
nected  to  a  computer  with  LabVIEW  software  (National  In¬ 
struments  Corporation,  Austin,  TX).  Post  processing  of  the  data  and 
the  nonlinear  analysis  were  done  with  custom  codes  run  using 


pressure  taps 


Fig.  6.  Custom,  fully-parallel  experimental  cathode  flow  field  with  pressure  taps  in  the 
headers. 


MATLAB  software  (MathWorks,  Natick,  MA).  Computational  pro¬ 
cessing  was  accelerated  by  writing  and  compiling  the  neighbor¬ 
hood  search  algorithms  in  C,  as  well  as  by  parallelizing  the  code. 

4.  Results 

To  characterize  the  overall  fuel  cell  performance,  we  measured 
galvanostatic  polarization  curves  at  the  experimental  conditions  for 
air  stoichiometric  flow  rates  of  ?  =  1.25, 1.5,  and  2.  The  VI  curves  are 
shown  in  Fig.  7.  At  low  current  densities,  there  is  a  reduction  in 
voltage  with  decreasing  ?  attributable  to  channel  flooding.  The 
presence  of  liquid  water  in  the  channels  at  low  current  densities  is 
typical  of  a  two-phase  slug  flow  regime,  and  decreases  the  effective 
active  area  of  the  fuel  cell  resulting  in  voltage  loss.  There  is  a  slight 
recovery  of  the  voltage  loss  for  ?  =  1.25  at  0.4  A  cm-2  as  the  amount 
of  liquid  water  present  in  the  channels  decreases.  At  high  current 
densities,  the  concentration  overpotential  becomes  evident.  The 
high-current  polarization  at  ?  =  1.25  may  also  have  effects  from  the 
presence  of  liquid  water  in  the  catalyst  and  gas  diffusion  layers, 
causing  increased  resistance  to  oxygen  diffusion. 

4.2.  Chaotic  dynamics 

We  conducted  galvanostatic  potentiometry  experiments  to 
identify  fuel  cell  operating  conditions  under  which  the  output 
voltage  signal  is  chaotic.  Due  to  the  degrees  of  freedom  provided  by 
the  parallel  air  delivery  microchannels,  we  expected  conditions 
that  lead  to  slug  and  plug  channel  flow  regimes  to  yield  the  most 
chaotic  dynamics.  Noting  the  scalings  ug  °c  ?i  and  u/  °c  i,  where  i  is 
the  current,  two-phase  flow  regime  maps  [4,5]  indicate  that  low- 
current,  low-?  conditions  are  most  chaotic,  consistent  with  our 
own  experimental  data. 

Fig.  8  shows  20  s  durations  of  the  voltage  signals  for  current 
densities  of  0.1  A  cm-2, 0.2  A  cm-2,  and  0.3  A  cm-2.  For  each  current 
density,  voltage  time  series  were  collected  for  air  stoichiometric 
ratios  of  ?  =  1.25, 1.5, 1.75,  and  2.  The  computed  nonlinear  statistics 
presented  later  are  based  on  500  s  durations  of  the  voltage  time 
series  after  the  equilibration  time.  At  a  constant  current  density,  the 
voltage  fluctuation  magnitude  decreases  as  ?  is  increased,  indi¬ 
cating  more  effective  removal  of  liquid  water  due  to  the  increased 


Fig.  7.  Galvanostatic  VI  curves  at  experimental  conditions  with  £  =  1.25  (bottom), 
1.5  (middle),  and  2  (top). 
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Fig.  8.  Voltage  time  series  at  £  =  1.25, 1.5, 1.75,  and  2  (bottom  to  top)  for  varying  current  densities  of  a)  0.1  A  cm  2,  b)  0.2  A  cm  2,  and  c)  0.3  A  cm  2.  Fluctuation  magnitude  decreases 
with  increasing  £,  and  fluctuation  frequency  increases  with  increasing  current. 


air  flow  rate  [8].  As  the  current  density  increases,  and  hence  the  air 
and  water  flow  rates,  so  does  the  voltage  fluctuation  frequency. 

For  the  same  experimental  conditions,  we  averaged  the  differ¬ 
ential  two-phase  cathode  pressure  drop,  AP20  (the  subscript  20  is 
used  to  denote  two-phase),  over  the  500-s  durations  used  for  the 
nonlinear  analysis  and  plotted  them  versus  ug  in  Fig.  9.  For  laminar, 
single-phase  flow  of  the  gas,  we  expect  a  linear  increase  in  APg  with 
ug  (the  subscript  g  denotes  a  gas-only  flow)  from  the  theory  of  flow  in 
pipes.  However,  for  two-phase  flow  the  pressure  drop  is  commonly 
modeled  as  A P2(p  =  A Pg02(p,  where  $20  is  the  two-phase  friction 
multiplier.  Due  to  the  nonlinearities  of  two-phase  flow,  <P 2 <p  is  an 
empirical  parameter  that  corrects  the  gas-only  pressure  drop  for  the 
presence  of  another  phase,  and  is  thus  dependent  upon  the  two- 
phase  flow  regime.  The  decrease  in  A P2^  with  increasing  ug  at 
0.1  A  cm”2  can  be  attributed  to  a  decrease  in  such  as  would  occur 
in  the  onset  of  the  transition  to  a  flow  regime  with  a  lower  pressure 
drop.  The  nonlinear  trend  at  0.2  A  cm”2  can  be  likewise  attributed  to 
a  decrease  in  <P20.  At  0.3  A  cm”2,  we  see  a  linear  trend  indicating  an 


Fig.  9.  Differential  two-phase  cathode  pressure  A P2(p  vs.  the  superficial  air  velocity  ug 
for  current  densities  of  0.1  A  cm  2,  0.2  A  cm”2,  and  0.3  A  cm”2,  each  with  £  =  1.25, 1.5, 
1.75,  and  2.  The  values  shown  are  averages  over  the  500-s  durations  used  for  the 
nonlinear  analysis,  and  the  error  bars  are  the  standard  deviations.  The  nonlinearities  of 
two-phase  flow  can  be  seen  in  the  0.1  A  cm  2  and  0.2  A  cm”2  differential  pressures, 
which  do  not  increase  linearly  with  ug  as  expected  for  single-phase  flow. 


approximately  constant  <f>2<p  and  hence  a  stable  two-phase  flow 
regime  like  annular  flow.  For  a  discussion  of  <P20,  see  Section  5.1. 

After  filtering  the  500-s  durations  of  the  voltage  signals  with  the 
nonlinear  noise  reduction  algorithm  from  Section  2.1,  we  computed 
correlation  sums  to  use  in  the  estimation  of  the  chaotic  invariants 
D2  and  I<2.  The  removal  of  sufficient  noise  to  reveal  power  law 
scaling  proved  to  be  sensitive  to  the  parameters  chosen  for  the 
noise  reduction,  which  require  some  a  priori  knowledge  of  the 
attractor  dimension  and  the  noise  level.  To  account  for  the  spurious 
effects  of  noise,  which  has  a  greater  impact  at  lower  scaling  ranges 
and  causes  overestimation  of  the  invariants  [29  ,  we  weighted  our 
estimates  linearly  by  the  median  value  of  the  (un-normalized) 
scaling  range  (for  a  discussion  of  the  weighting,  see  Section  5.1). 

Fig.  10  plots  the  weighted  D2  and  K2  estimates  as  a  function  of 
current  density  and  J,  where  the  error  bars  cover  the  variability  in 
C(r,m)  over  the  scaling  range.  The  conditions  for  which  no  estimates 
are  shown  are  those  that  did  not  yield  convincing  power  law  scaling 
in  the  correlation  sums.  This  estimate  uncertainty  could  result  from 
an  insufficient  amount  of  data,  the  presence  of  noise,  or  the  lack  of 
chaotic  behavior. 

The  positive  entropy  JC2  is  a  sufficient  condition  to  indicate 
chaotic  dynamics  [24],  and  the  fractal  dimension  D2  further  sup¬ 
ports  this  conclusion  while  providing  us  with  the  embedding  di¬ 
mensions  required  to  compute  accurate  nonlinear  statistics.  We 
observed  that  the  weighted  D2  and  K2  estimates  for  0.2  and 
0.3  A  cm”2  fall  on  the  same  curve,  while  those  for  0.1  A  cm”2 
remain  slightly  higher,  at  least  for  f  =  1.5  and  1.75.  The  estimates 
show  a  reduction  in  the  active  degrees  of  freedom  and  the  entropy 
for  increasing  f ,  indicating  that  the  dynamics  are  more  stable  with 
higher  air  flow  rates  consistent  with  more  effective  cathode  water 
removal.  The  small  reduction  in  D2  and  I<2  from  0.1  A  cm”2  to 
0.2  A  cm”2  indicates  that  for  the  conditions  investigated,  the  added 
stability  from  the  increase  in  the  air  flow  rate  corresponding  to  the 
increase  in  current  density  (for  constant  f )  is  mostly  offset  by  the 
increase  in  water  production.  This  ability  to  quantify  the  effect  of 
operating  conditions  on  dynamical  stability  associated  with  cath¬ 
ode  water  accumulation  can  be  used  to  stabilize  PEFC  operation. 

Unfortunately,  nonlinear  noise  reduction  and  computing  the 
correlation  sums  for  invariant  estimation  are  computationally 
intensive.  For  instance,  on  a  performance  workstation,  D2  and  K2 
estimates  require  at  least  10  times  the  real-time  length  of  the  signal 
using  code  that  has  been  optimized  for  computational  speed.  Thus, 
the  computational  time  does  not  make  them  amenable  to  real¬ 
time,  embedded  control.  Even  so,  the  D2  and  K2  parameters  pro¬ 
vide  a  richness  of  results  that  is  valuable  for  post-experimental  data 
analysis  and  are  well-suited  to  that  usage. 
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Fig.  10.  a)  Correlation  dimension  and  b)  Kolmogorov  entropy  estimates  for  varying  £  at 
constant  current  densities  of  0.1  A  cm  2,  0.2  A  cnrr2,  and  0.3  A  cm  2.  Estimates  were 
weighted  by  the  scaling  range,  with  unweighted  D2  estimates  ranging  from  about  7  to 
14.  The  error  bars  cover  the  variability  in  C(r,m)  over  the  scaling  range. 

4.2.  Lyapunov  exponents 

For  embedded  PEFC  control  with  nonlinear  diagnostics,  it  is  vital 
that  the  nonlinear  signal  analysis  is  not  computationally  burden¬ 
some  in  order  to  allow  low-cost  electronics.  As  an  alternative 
approach  to  D2  and  I< 2  computation,  the  unfiltered  voltage  signals 
were  converted  to  one-dimensional  return  maps  to  evaluate 
whether  the  reduced  order  map  of  the  raw  voltage  signal  main¬ 
tained  some  of  the  chaotic  properties  of  the  continuous  time  state- 
space  attractor  of  the  noise-reduced  voltage  signal.  Voltage  peaks 
were  taken  to  be  at  least  a  minimum  percentage  of  the  standard 
deviation  of  the  data  above  the  previous  valley  (the  tolerance),  and 
the  Lyapunov  exponents  of  the  return  map  were  estimated  from 
the  mean  divergence  of  the  50  nearest  peaks.  For  the  values  re¬ 
ported  in  Fig.  11,  we  averaged  the  LEs  computed  from  peaks  over  a 
representative  tolerance  range  of  40%-60%,  and  the  error  bars  are 
the  standard  deviations  of  the  LEs  from  the  same  range. 

The  Lyapunov  exponents  of  the  voltage  return  map  in  Fig.  11 
show  the  same  trend  as  the  rigorous  D2  and  I<2  estimates  of  the 
time-delay  embedded  state  space.  Lyapunov  exponents  decrease 
with  increasing  ?,  and  the  curve  for  0.1  A  cm-2  is  different  than  that 
for  0.2  and  0.3  A  cm-2.  Additionally,  the  LE  analysis  provides  esti¬ 
mates  for  the  conditions  at  which  D2  and  K 2  estimation  was  not 
possible  due  to  uncertainty.  While  the  Lyapunov  exponents  of  the 
voltage  return  map  are  reduced-order  projections  of  the  Lyapunov 
exponents  from  the  full  state  space,  the  comparison  with  D2  and  I< 2 
shows  that  they  capture  the  relevant  stability  characteristics  of  the 
dynamics  and  thus  are  an  accurate  indication  of  PEFC  stability.  The 
ease  of  calculation  of  the  map  Lyapunov  exponents  compared  to  D2 
and  J<2  make  them  excellent  candidates  for  real  time  stability 
control,  discussed  further  in  Section  5.2. 
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Fig.  11.  Lyapunov  exponent  estimates  of  the  voltage  return  map  for  varying  q  at 
constant  current  densities  of  0.1  A  cnrr2,  0.2  A  cnrr2,  and  0.3  A  cm-2.  The  estimates  are 
averaged  from  the  mean  value  of  the  50  closest  peaks  for  varying  peak  tolerances,  and 
the  error  bars  are  the  standard  deviations. 

The  large  standard  deviation  of  the  local  LE  estimates,  shown  by 
the  error  bars  in  Fig.  11,  makes  it  necessary  to  include  a  sufficiently 
large  number  of  voltage  peaks  in  order  to  converge  to  the  correct  LE 
value.  To  determine  how  many  voltage  peaks  are  required,  we 
computed  LE  estimates  for  0.1  A  cm-2  over  varying  time  windows 
of  data.  Fig.  12  shows  the  convergence  of  the  LE  estimates  for 
increasing  time  window  sizes,  where  the  shaded  error  is  the 
standard  deviation  of  the  LE  estimates  of  the  previous  10  time 
windows.  It  can  be  seen  that  for  time  windows  >200  s  the 
?  =  1.25  LE  estimate  begins  to  be  distinguished  from  the  ?  =  1.5  LE 
estimate,  with  convergence  to  the  appropriate  value  being  reached 
around  300  s.  For  the  less  extreme  cases  of  ?  =  1.5, 1.75,  and  2,  the 
LE  estimates  converge  much  more  quickly  and  are  distinct  for  all 
time  windows  shown. 

4.3.  Hurst  exponents 

To  gain  more  insight  into  how  the  fuel  cell  dynamics  change 
with  the  investigated  operating  conditions,  we  used  diffusion 
analysis  of  350  s  of  the  raw  voltage  signals  to  estimate  the  Hurst 
exponents.  Fig.  13  shows  that  for  all  the  conditions,  the  shortest 
time  windows  yield  local  power  law  exponents  >0.9,  indicative  of  a 
strong,  positive  autocorrelation.  This  can  be  attributed  to  the 
redundancy  in  the  highly-sampled  data  and  the  strong  deter¬ 
minism  over  short  times.  As  the  time  window  increases,  the  local 
power  law  exponents  decrease  until  time  windows  large  enough  to 
capture  constant  power  law  scaling  in  the  voltage  diffusion  are 
used.  The  power  law  scaling  is  limited  at  the  longest  time  windows 
by  curvature  and  kinks  in  the  diffusion  analysis,  which  may  be 
spurious  artifacts  resulting  from  finite  data  [30]. 

The  onset  of  power  law  scaling,  seen  as  linearity  in  the  log-log 
plot,  as  well  as  the  Hurst  exponents  estimated  in  the  scaling  range, 
are  dependent  upon  operating  conditions.  The  time  window  length 
where  power  law  scaling  begins  is  longest  for  low-current,  low-? 
conditions,  and  decreases  with  increasing  current  and  ?  as  the 
dynamical  timescales  decrease  due  to  increasing  air  and  water  flow 
rates.  The  Hurst  exponent  estimates  at  0.1  A  cm-2  and  low  ?  are 
close  to  H  =  0.5,  indicating  a  lack  of  correlation  due  to  the  high 
entropy  and  instability  in  the  chaotic  dynamics  from  a  significant 
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Fig.  12.  Voltage  return  map  Lyapunov  exponent  estimates  at  0.1  A  cm  2  vs.  the  time 
window  used  for  the  estimate.  The  extreme  low  current/low  £  case  shows  that  be¬ 
tween  200  and  300  s  of  data  are  needed  to  get  a  conservative  50  peaks  with  neighbors 
that  accurately  represent  the  dynamics.  The  shaded  error  is  the  standard  deviation  of 
the  LE  estimates  of  the  previous  10  time  windows. 

amount  of  water  in  the  cathode.  As  the  current  and  ?  increase,  the 
increasing  value  of  H  shows  longer  autocorrelations  as  the  dy¬ 
namics  gain  memory  from  increasing  stability. 

In  the  diffusion  analysis,  the  most  chaotic  conditions  have  weak 
autocorrelations  and  appear  stochastic  at  sufficiently  long  time 
windows,  and  as  the  chaos  subsides  the  strength  of  the  autocor¬ 
relation  increases.  Like  we  observed  in  the  D2  and  K 2  analysis, 
0.2  A  cm-2  and  0.3  A  cnrT2  show  similar  behavior  exhibited  by 
similar  power  law  scaling,  with  the  exception  of  0.3  A  cm-2,  ?  =  2. 
At  this  condition,  no  chaotic  invariants  were  estimated  and  the 
diffusion  analysis  reveals  a  higher  H,  indicating  a  change  in 
dynamical  behavior  from  the  transition  to  a  more  stable  operating 
regime.  The  diffusion  analysis  provides  insight  into  the  timescales 
relevant  to  PEFC  dynamical  stability  as  well  as  an  additional  mea¬ 
sure  of  the  strength  of  the  stability  for  the  operating  conditions 


Fig.  13.  Diffusion  analysis  for  varying  £  at  constant  current  densities  of  0.1  A  cm  2, 
0.2  A  cm-2,  and  0.3  A  cnrr2.  The  curves  have  been  vertically  separated  for  clarity. 


investigated.  Both  of  these  results  are  important  for  controller 
design,  discussed  further  in  Section  5.2. 

5.  Discussion 

The  presence  of  high-dimensional  chaos  in  low-current,  low-? 
PEFC  operation  indicates  the  complexity  and  stability  of  the  dy¬ 
namics  governing  cathode  water  accumulation  and  removal, 
including  two-phase  flow  in  parallel  microchannels.  The  estimation 
of  invariants  D2  and  K 2  is  a  rigorous  way  of  quantifying 
experimentally-obtained  chaotic  PEFC  dynamics  and  the  related 
stability.  However,  there  are  some  drawbacks,  namely  the  compu¬ 
tational  power  and  tuning  of  parameters  required  to  obtain  accu¬ 
rate  estimates,  that  make  D2  and  I< 2  difficult  to  compute  as  part  of  a 
quick,  robust  PEFC  water  management  control  scheme.  The  D2  and 
K 2  results  presented  in  this  paper  inform  the  development  of  a 
reduced-order  nonlinear  statistic,  the  Lyapunov  exponent  of  the 
voltage  return  map,  that  provides  a  first-order  approximation  of 
the  robust  invariants  while  requiring  minimal  effort  to  compute.  In 
this  section,  we  discuss  the  D2  and  K 2  estimation  procedure  and  its 
drawbacks,  the  chaotic  dynamics  and  their  relation  to  cathode 
water  content,  and  the  Lyapunov  and  Hurst  exponents  and  their 
applications  to  active  PEFC  control. 

5.2.  Chaotic  behavior:  D2  and  I<2  estimates 

The  estimation  of  chaotic  invariants  from  an  experimental  time 
series  requires  a  sufficiently  long  period  of  data  to  populate  the 
attractor  to  the  point  where  power  law  scaling  is  revealed.  As  the 
number  of  degrees  of  freedom  increases,  the  presence  of  positive 
Lyapunov  exponents  becomes  more  apparent,  making  it  necessary 
to  obtain  exponentially  increasing  amounts  of  data  to  make  accu¬ 
rate  estimates.  To  obtain  a  scaling  range  of  at  least  a  decade,  it  is 
estimated  that  more  than  10D2/2  points  are  required  [31  .  The  un¬ 
weighted  D2  estimates  for  the  PEFC  dynamics  studied  ranged  from 
about  7  to  14  and  required  large  data  sets.  Without  using  the 
nonlinear  noise  reduction,  our  correlation  sums  showed  no  power 
law  scaling.  Even  with  the  noise  reduction,  our  scaling  ranges  were 
very  small  (less  than  a  decade)  and  noise  still  had  a  significant 
presence  in  the  voltage  signal.  However,  this  is  expected  for  such 
high-dimensional  chaos  with  a  limited  amount  of  data  [27].  In 
addition,  the  definite  power  law  scaling  in  the  correlation  sums 
(Fig.  4)  demonstrates  the  reliability  of  our  estimates.  The  observed 
D2  trend  in  Fig.  10  shows  an  increasing  D2  with  decreasing  current 
and  ?,  making  it  likely  that  the  low-current,  low-?  conditions  with 
no  estimates  (e.g.  0.1  A  cm-2,  ?  =  1.25)  did  not  have  enough  data 
points  to  reveal  the  power  law  scaling.  For  the  missing  estimate  at 
the  highest  current  and  ?  tested  (0.3  A  cm-2  and  2,  respectively),  it 
is  possible  that  the  dynamical  behavior  changed  in  the  transition  to 
a  stable  channel  flow  regime,  as  it  is  unlikely  that  chaotic  behavior 
ceased  altogether  at  this  condition.  Such  a  change,  which  can  be 
seen  in  the  Hurst  exponent  estimates  in  Fig.  13,  could  alter  the 
dimensionality  and  scaling  regions  of  the  attractor,  rendering  our 
nonlinear  noise  reduction  parameters  ineffective. 

As  discussed  earlier,  the  presence  of  noise  in  the  signal  has  been 
shown  to  cause  overestimation  of  D 2  [29  .  To  better  understand  the 
effects  of  noise  on  chaotic  dynamics,  we  computed  correlation 
sums  from  the  well-characterized  chaotic  Lorenz  system  [20]  with 
superposed  white  Gaussian  measurement  and  dynamical  noise. 
The  correlation  sums  of  the  noisy  Lorenz  system  also  demonstrated 
increasing  overestimation  of  D2  at  decreasing  scaling  ranges. 
Schreiber  has  shown  that  correlation  dimensions  estimated  from 
correlation  sums  with  measurement  noise  follow 
D2(r,m  +  1 )  =  D2(r,m)  +  D2>noise(r)  [32].  For  Gaussian  noise,  it  can  be 
shown  that  D2,Gauss(r)  =  (r  exp(-r2/4ff2))/(ff^erf(r/2<T)),  a 
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decreasing  function  of  r  (here  a  is  the  noise  amplitude).  So,  to  ac¬ 
count  for  noise  in  our  D2  and  I< 2  estimates,  we  weighted  the  esti¬ 
mates  linearly  by  the  median  value  of  the  scaling  range,  with  higher 
weight  given  to  larger  scaling  ranges.  While  the  linear  weighting  is 
not  exact  for  Gaussian  noise,  it  serves  as  a  first  order  approxima¬ 
tion.  This  weighting  collapsed  the  0.2  A  cnrT2  and  0.3  A  cm-2  D2  and 
I<2  estimates  onto  a  single  curve  while  the  0.1  A  cm”2  estimates 
remained  higher  (see  Fig.  10).  The  unweighted  Lyapunov  exponent 
estimates  of  the  unfiltered  voltage  return  map  show  a  similar  trend, 
Fig.  11,  as  do  the  Hurst  exponent  estimates,  Fig.  13,  indicating  that 
all  of  the  statistics  are  describing  similar  nonlinear  phenomena. 

The  collapse  of  the  0.2  A  cm”2  and  0.3  A  cm”2  chaotic  invariants 
imply  that  for  this  experimental  setup,  chaotic  behavior  is  more 
sensitive  to  changes  in  air  stoichiometric  ratio  (or  superficial  gas 
velocity  ug)  than  to  changes  in  current  (or  superficial  liquid  velocity 
ui)  beyond  a  threshold  current  density  between  0.1  A  cm”2  and 
0.2  A  cm”2.  This  behavior  is  indicative  of  a  bifurcation  in  the  dy¬ 
namics,  a  nonlinear  phenomena  whereby  the  attractor  topology 
will  suddenly  change,  sometimes  dramatically,  over  a  smooth 
change  in  the  system  parameters.  In  a  PEFC,  such  a  bifurcation 
could  result  from  a  change  in  the  two-phase  flow  regime. 

As  an  indicator  of  the  two-phase  flow  regime,  we  calculated  the 
two-phase  friction  factor  &2 0  such  that  AP20  =  kPg$2(p  where  A Pg 
was  recorded  for  the  dry  cathode  flow  field  with  no  current.  Fig.  14 
shows  the  LE  estimates  versus  <P2(p.  The  high  &2(p  for  the  0.1  A  cm”2 
conditions  indicates  significant  two-phase  flow  effects  arising  from 
slug  and  plug  flow  when  ug  and  u/  are  small.  As  the  current  density 
and  ?  increase,  $20  approaches  unity  and  the  two-phase  effects 
subside  in  the  transition  to  a  regime  with  less  interaction  between 
phases  like  the  annular/film  flow  regime.  The  0.1  A  cm”2  conditions 
result  in  a  u/  value  of  6  x  10”5  m  s”1  and  ug  values  from  0.3  to 
0.5  m  s”1  that  are  well  inside  the  slug  flow  regime  outlined  by 
Hussaini  and  Wang  [4]  and  exhibit  the  corresponding  high  two- 
phase  friction  multipliers  expected  in  this  flow  regime.  The  high- 
?  conditions  for  the  0.2  and  0.3  A  cm”2  conditions  result  in  super¬ 
ficial  velocities  approaching  the  transition  to  film  flow,  and  as  ex¬ 
pected  the  two-phase  friction  multiplier  approaches  one.  At 
0.3  A  cm”2  with  f  =  1.75  and  2,  the  superficial  velocities  u/  and  ug 
are  within  the  film  flow  regime  bounds  measured  by  Hussaini  and 
Wang  and  result  in  two-phase  friction  multipliers  of  about  one. 
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Fig.  14.  Two-phase  flow  regime  indication  based  on  The  0.1  A  cm  2  conditions 
appear  to  produce  a  different  flow  regime  as  shown  by  the  high  <t>2<p-  Convergence  to  a 
common  LE  estimate  may  be  seen  as  $2 <p  approaches  1. 


We  see  a  smooth  decrease  in  <P2(p  during  the  transition  rather 
than  an  abrupt  change  because  the  product  water  in  a  PEFC  cathode 
flow  field  is  introduced  along  the  length  of  the  channel.  The  water 
production  along  the  channel  causes  zi/  to  vary  along  the  length  of 
the  channel.  Thus,  for  outlet  values  of  u/  and  ug  near  the  slug  flow  to 
film  flow  transition,  the  flow  regime  can  also  vary  along  the  length 
of  the  channel,  yielding  a  convolution  of  the  distinct  flow  regime 
dynamics  in  the  observed,  overall  dynamics.  The  relative  amount  of 
each  regime  present  in  the  flow  field  depends  on  the  proximity  of 
the  flow  conditions  to  that  of  the  regime  transition.  The  experi¬ 
mental  ug  and  ui  values  for  these  observations  are  in  accordance 
with  the  flow  regime  map  and  transitions  reported  in  Ref.  [4]. 

As  the  two-phase  flow  effects  subside  and  $2 0  approaches  unity, 
the  LE  estimates  approach  a  value  around  3,  implying  the  presence 
of  a  LE  floor  due  to  dynamics  unrelated  to  the  two-phase  flow.  The 
magnitude  of  departure  from  the  LE  floor  with  decreasing  f  has 
some  proportionality  to  02 0,  but  we  expect  that  there  are  also 
contributions  from  concentration  overpotential.  With  increasing 
current  density  these  contributions  become  more  important,  and 
hence  we  see  minimal  difference  between  the  LEs  for  0.2  and 
0.3  A  cm”2  even  though  there  are  differences  in  $20  at  these  con¬ 
ditions.  The  relationship  between  the  Lyapunov  exponents  and  &2 0 
in  Fig.  14  shows  that  the  Lyapunov  exponents  may  be  calculated  as 
part  of  a  control  scheme  to  provide  information  about  the  cathode 
water  content  as  it  relates  to  the  two-phase  flow  regime  and  PEFC 
dynamical  stability,  as  discussed  below. 

5.2.  Control  applications:  Lyapunov  and  Hurst  exponent  estimates 

For  efficient  and  effective  PEFC  water  management,  there  are 
potentially  great  benefits  from  having  nonlinear  diagnostics  indi¬ 
cate  the  health  of  operating  PEFCs.  Control  systems  based  on  linear 
behavior  like  fluctuation  magnitude  may  take  corrective  action 
when  none  is  needed.  For  example,  the  standard  deviation  of  the 
0.1  A  cm”2,  ?  =  1.75  voltage  signal  is  about  half  that  of  the 
0.3  A  cm”2,  f  =  1.5  voltage  signal  (Fig.  8),  while  the  D2  estimate  is 
greater  and  the  K 2  estimate  about  the  same  (Fig.  10).  Thus  it  may 
prove  beneficial  to  base  corrective  action  on  chaotic  properties  like 
D2  and  K 2.  However,  the  estimation  procedure  for  these  invariants 
requires  a  large  number  of  data  points  and  is  computationally 
expensive,  as  mentioned  earlier,  making  them  unlikely  candidates 
for  real-time  control. 

Instead,  we  propose  using  the  Lyapunov  exponents  of  the 
voltage  return  map  to  indicate  chaotic  instability  resulting  from 
problematic  cathode  water  accumulation.  Elevation  of  the  Lyapu¬ 
nov  exponents  from  a  known  stable  condition,  which  may  vary 
between  fuel  cells,  shows  that  the  two-phase  flow  is  approaching 
an  unstable  regime  and  corrective  action  should  be  taken.  Fig.  14 
shows  that  the  LE  estimates  converge  to  a  floor  value  as  the  two- 
phase  flow  becomes  stable,  so  a  stable,  reference  LE  value  in  a 
control  scheme  can  be  selected  from  such  a  plot.  Use  of  the  return 
map  reduces  the  order  of  the  analysis  from  as  many  as  14  di¬ 
mensions  to  1  dimension.  While  such  a  large  reduction  of  order 
may  eliminate  many  dynamical  properties,  we  have  shown  that  the 
Lyapunov  exponent  estimates  of  the  unfiltered  voltage  return  map 
follow  the  same  trend  as  the  weighted  D2  and  K 2  estimates  of  the 
noise-reduced  embedded  attractor.  The  average  inter-peak  in¬ 
tervals  used  for  the  return  map  range  from  0.1  s  (high  current  and  f ) 
to  1.7  s  (low  current  and  f ).  Hence,  the  number  of  time  series  points 
used  for  LE  estimation  is  one  to  two  orders  of  magnitude  less  than 
that  used  for  D2  and  K 2  estimation.  Using  a  conservative  0(log  n) 
neighborhood  search  complexity  and  taking  into  account  the  total 
number  of  neighborhood  searches  required,  the  LE  estimation 
procedure  is  two  to  three  orders  of  magnitude  faster  than  the  D2 
and  I<2  estimation  procedure  from  the  reduction  in  data  alone.  The 
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Lyapunov  exponents  require  only  seconds  of  serial  computation  for 
estimation,  while  D2  and  K 2  estimates  require  hours  of  parallel 
processing. 

Fig.  1  shows  a  300  s  time  lag  between  the  onset  of  chaotic 
behavior  and  complete  voltage  shutdown  from  flooding.  In  an 
automotive  stack  with  constantly  changing  parameters,  300  s 
might  be  too  large  a  window  to  provide  effective  cathode  flooding 
mitigation.  It  can  be  seen  in  Fig.  12  that  most  conditions  can  be 
distinguished  using  the  LE  analysis  in  as  little  as  30  s,  the  shortest 
time  window  shown.  Furthermore,  this  analysis  uses  a  conservative 
50  peaks  for  the  LE  estimate.  By  implementing  a  minimum  sepa¬ 
ration  for  peaks  to  be  considered  neighbors,  it  is  possible  to  achieve 
an  accurate  estimate  with  fewer  peaks.  When  coupled  with 
knowledge  of  the  distributions  of  local  LEs,  changes  in  dynamical 
behavior  could  be  detected  in  seconds  from  a  few  peaks  indicative 
of  the  dynamics. 

Another  candidate  metric  for  chaotic  instability  detection  is  the 
Hurst  exponent  estimate  from  the  diffusion  analysis  of  the  voltage 
signal.  The  diffusion  analysis  in  Fig.  13  shows  some  interesting 
properties  of  high-dimensional  chaos,  as  well  as  relevant  time- 
scales  of  the  dynamics  for  our  PEFC  under  the  operating  conditions 
investigated.  In  the  high-dimensional  limit,  deterministic  chaotic 
dynamics  may  appear  stochastic  [33  .  This  is  due  to  the  positive 
Lyapunov  exponent  present  in  deterministic  chaos.  Small  errors  are 
compounded  exponentially  quickly,  and  when  many  degrees  of 
freedom  are  present  it  will  cause  the  dynamics  at  large  enough  time 
windows  to  appear  stochastic.  From  the  diffusion  analysis  of  our 
PEFC  system,  we  can  see  that  the  departure  from  short-time,  bal¬ 
listic  determinism  occurs  after  approximately  1  s.  Interestingly,  this 
is  also  the  timescale  of  the  first  minimum  in  the  mutual  informa¬ 
tion  of  the  signal,  which  we  used  as  our  embedding  time  delay, 
typically  around  0.5  s. 

The  results  of  the  voltage  diffusion  analysis,  specifically  the 
power  law  scaling  transitions  between  operating  conditions  and 
time  windows,  can  be  used  for  PEFC  water  management.  Dynamics 
are  monitored  for  stochastic  behavior  at  time  intervals  >1  s, 
through  Hurst  exponent  estimation  or  through  the  use  of  nonlinear 
prediction.  As  the  dynamics  appear  more  stochastic,  the  Hurst 
exponent  will  approach  0.5  and  the  predictability  will  worsen. 
When  a  predefined  threshold  is  crossed,  corrective  action  could  be 
taken. 

In  future  work,  the  effects  of  PEFC  operating  conditions,  such  as 
inlet  relative  humidity  and  temperature  as  well  as  flow  field  design 
and  wettability  should  be  explored  to  elucidate  their  effects  on 
PEFC  stability. 

6.  Conclusion 

In  this  research,  we  have  detected  the  presence  of  high¬ 
dimensional,  chaotic  dynamics  in  PEFC  operation  at  low-current, 
low-?  conditions  typical  of  startup  and  idling.  The  chaotic 
behavior  increases  with  decreasing  current  and  decreasing  f ,  as  the 
channel  ug  and  u/  move  towards  those  of  the  two-phase  slug  flow 
regime.  By  weighting  our  chaotic  invariant  estimates  of  dimension 
and  entropy  by  the  median  of  the  state-space  power  law  scaling 
range  to  account  for  noise,  we  detected  bifurcation-like  dynamics 
consistent  with  a  two-phase  flow  regime  transition. 

We  identified  a  reduced-order  metric  indicative  of  the  chaotic 
behavior.  The  Lyapunov  exponents  of  the  one-dimensional  voltage 
return  map  capture  the  chaotic  behavior  of  the  high-dimensional 
reconstructed  state  space  as  quantified  by  the  correlation  dimen¬ 
sion  and  Kolmogorov  entropy.  The  inclusion  of  >300  s  of  data 
provides  an  accurate  LE  estimate  for  all  conditions  investigated.  For 
use  in  a  PEFC  water  management  scheme,  effective  LE  estimation 
from  the  voltage  return  map  is  achievable  in  seconds. 


Using  a  diffusion  analysis  to  estimate  the  Hurst  exponents  of  the 
voltage  signal,  we  have  seen  how  the  high-dimensional  chaos  af¬ 
fects  the  long-term  memory  of  the  PEFC  dynamics.  More  chaotic 
dynamics  have  weaker  autocorrelations  at  long  times.  Detection  of 
stochastic  behavior  at  a  time  lag  >1  s  can  indicate  the  level  of 
chaotic  behavior  and  be  applied  in  PEFC  control. 
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